Ion cyclotron instability

note: Using 1-d full PIC code written in FORTRAN (other files).

parameters

  • $m_i/m_e=25$
  • $\omega_{pe}/\Omega_{ce}=5$
  • $n_b/n_0=0.2$
  • $v_d/v_A=3.0$
  • $\beta_e=0.5$
  • $\beta_i=0.05$
  • $ppc=200$

Overview

In [1]:
%config InlineBackend.figure_format='retina'
import numpy as np
import matplotlib.pyplot as plt
import math as mt
from scipy.fftpack import fft, ifft,fft2,ifft2,fftshift,ifftshift
from IPython import display
import matplotlib.animation as animation
In [2]:
plt.rcParams['animation.html'] = 'html5'
In [3]:
bs=1024
nb=0.2

#fig,(ax1,ax2)=plt.subplots(1,3,figsize=(10,5))
fig=plt.figure(figsize=(12,4))
ax1=plt.subplot2grid((1,3),(0,0),colspan=2)
ax2=plt.subplot2grid((1,3),(0,2))

def update_anim(i):
    it=i*16
    f1=open('xp1_%05d' %(it),'rb') 
    f2=open('xp2_%05d' %(it),'rb')
    f3=open('xp3_%05d' %(it),'rb')
    xp1=np.fromfile(f1,dtype='float64').reshape(-1,4)
    xp2=np.fromfile(f2,dtype='float64').reshape(-1,4)
    xp3=np.fromfile(f3,dtype='float64').reshape(-1,4)
    
    for ax in (ax1,ax2): 
        ax.clear()

    ax1.plot(xp1[:,0],xp1[:,1],'C0,')
    ax1.plot(xp2[:,0],xp2[:,1],'C1,')
    ax1.plot(xp3[:,0],xp3[:,1],'C2,')
    ax1.set_xlabel(r'$x/(c/\omega_{pe})$')
    ax1.set_ylabel(r'$v_x/c$')
    ax1.set_ylim(-0.5,0.5)

    ax2.hist(xp1[:,1],bins=bs,range=(-0.25,0.25),histtype='step',color='b',weights=(1-nb)*np.ones(len(xp1)))
    ax2.hist(xp2[:,1],bins=bs,range=(-0.25,0.25),histtype='step',color='orange')
    ax2.hist(xp3[:,1],bins=bs,range=(-0.25,0.25),histtype='step',color='g',weights=nb*np.ones(len(xp1)))
    ax2.set_ylabel(r'$v_x/c$')
    ax2.set_xlabel(r'$f(v_x)$')
    ax2.set_ylim(0,6000)
    
anim=animation.FuncAnimation(fig,update_anim,frames=64)    
plt.close()
In [4]:
anim
Out[4]:

Time-space evolution of electromagnetic fields

In [5]:
%%time
dx=0.1
dt=0.1
isav=32
nt=4096
nx=1024
wce=0.2
rm=25
for it in range(1,3):
    f1=open('ex%03d.hst' %(it),'rb');f2=open('ey%03d.hst' %(it),'rb');f3=open('ez%03d.hst' %(it),'rb')
    f4=open('by%03d.hst' %(it),'rb');f5=open('bz%03d.hst' %(it),'rb')
    ex=np.fromfile(f1,dtype='float64').reshape(-1,nx)
    ey=np.fromfile(f2,dtype='float64').reshape(-1,nx)
    ez=np.fromfile(f3,dtype='float64').reshape(-1,nx)
    by=np.fromfile(f4,dtype='float64').reshape(-1,nx)
    bz=np.fromfile(f5,dtype='float64').reshape(-1,nx)
    if it==1: 
        exsave=ex;eysave=ey;ezsave=ez;bysave=by;bzsave=bz
    if it>1:
        exsave=np.r_[exsave,ex]
        eysave=np.r_[eysave,ey]
        ezsave=np.r_[ezsave,ez]
        bysave=np.r_[bysave,by]
        bzsave=np.r_[bzsave,bz]

CPU times: user 48.7 ms, sys: 208 ms, total: 257 ms
Wall time: 2.53 s
In [6]:
%%time
xmax=nx*dx
tmax=nt*isav*dt/rm*wce
fig=plt.figure(figsize=(20,10))
ax1 = fig.add_subplot(231)
ax2 = fig.add_subplot(232, sharex=ax1, sharey=ax1)
ax3 = fig.add_subplot(233, sharex=ax1, sharey=ax1)
ax4 = fig.add_subplot(235, sharex=ax1, sharey=ax1)
ax5 = fig.add_subplot(236, sharex=ax1, sharey=ax1)
im1=ax1.imshow(exsave,aspect='auto',origin='lower',cmap='terrain_r',extent=[0,xmax,0,tmax])
im2=ax2.imshow(eysave,aspect='auto',origin='lower',cmap='bwr',extent=[0,xmax,0,tmax],clim=(-0.01,0.01))
im3=ax3.imshow(ezsave,aspect='auto',origin='lower',cmap='bwr',extent=[0,xmax,0,tmax],clim=(-0.01,0.01))
im4=ax4.imshow(bysave,aspect='auto',origin='lower',cmap='bwr',extent=[0,xmax,0,tmax],clim=(-0.1,0.1))
im5=ax5.imshow(bzsave,aspect='auto',origin='lower',cmap='bwr',extent=[0,xmax,0,tmax],clim=(-0.1,0.1))
ax1.set_xlabel(r'$x/(c/\omega_{e})$');ax4.set_xlabel(r'$x/(c/\omega_{e})$');ax5.set_xlabel(r'$x/(c/\omega_{e})$')
ax1.set_ylabel(r'$t\Omega_{ci}$');ax4.set_ylabel(r'$t\Omega_{ci}$')
ax1.set_title(r'$E_x$');ax2.set_title(r'$E_y$');ax3.set_title(r'$E_z$')
ax4.set_title(r'$B_y$');ax5.set_title(r'$B_z$');
plt.colorbar(im1,ax=ax1)
plt.colorbar(im2,ax=ax2)
plt.colorbar(im3,ax=ax3)
plt.colorbar(im4,ax=ax4)
plt.colorbar(im5,ax=ax5)
plt.show()

CPU times: user 5 s, sys: 308 ms, total: 5.31 s
Wall time: 5.32 s

Time evolution of electromagnetic fields

In [7]:
nt=4096
enex=np.zeros(nt)
eney=np.zeros(nt)
enez=np.zeros(nt)
enby=np.zeros(nt)
enbz=np.zeros(nt)
for it in range(nt):
    enex[it]=np.average(exsave[it,:]**2)
    eney[it]=np.average(eysave[it,:]**2)
    enez[it]=np.average(ezsave[it,:]**2)
    enby[it]=np.average(bysave[it,:]**2)
    enbz[it]=np.average(bzsave[it,:]**2)

tt=np.linspace(0,nt*dt/rm*wce,nt)
plt.plot(tt,enex)
plt.plot(tt,eney)
plt.plot(tt,enez)
plt.plot(tt,enby)
plt.plot(tt,enbz)
plt.xlabel(r'$t\Omega_{ci}$')
plt.ylabel(r'$E_x^2$, $E_y^2$, $E_z^2$, $B_y^2$, $B_z^2$')
plt.yscale('log')
plt.ylim(1e-7,1e-2)
plt.show()

Time-space evolution of fluid values

In [8]:
%%time
bs=1024
nt=1024
den1save=np.zeros((nt,bs))
den2save=np.zeros((nt,bs))
den3save=np.zeros((nt,bs))
avl1save=np.zeros((nt,bs))
avl2save=np.zeros((nt,bs))
avl3save=np.zeros((nt,bs))
tmp1save=np.zeros((nt,bs))
tmp2save=np.zeros((nt,bs))
tmp3save=np.zeros((nt,bs))
for it in range(nt):
    f1=open('xp1_%05d' %(it),'rb') 
    f2=open('xp2_%05d' %(it),'rb')
    f3=open('xp3_%05d' %(it),'rb')
    xp1=np.fromfile(f1,dtype='float64').reshape(-1,4)
    xp2=np.fromfile(f2,dtype='float64').reshape(-1,4)
    xp3=np.fromfile(f3,dtype='float64').reshape(-1,4)

    (den1,bns)=np.histogram(xp1[:,0],bins=bs)
    (den2,bns)=np.histogram(xp2[:,0],bins=bs)
    (den3,bns)=np.histogram(xp3[:,0],bins=bs)
    (avl1,bns)=np.histogram(xp1[:,0],bins=bs,weights=xp1[:,1])
    (avl2,bns)=np.histogram(xp2[:,0],bins=bs,weights=xp2[:,1])
    (avl3,bns)=np.histogram(xp3[:,0],bins=bs,weights=xp3[:,1])
    tmp1=0.5*((xp1[:,1]-avl1[np.int_(xp1[:,0])]/den1[np.int_(xp1[:,0])])**2+xp1[:,2]**2+xp1[:,3]**2)
    tmp2=0.5*((xp2[:,1]-avl2[np.int_(xp2[:,0])]/den2[np.int_(xp2[:,0])])**2+xp2[:,2]**2+xp2[:,3]**2)
    tmp3=0.5*((xp3[:,1]-avl3[np.int_(xp3[:,0])]/den3[np.int_(xp3[:,0])])**2+xp3[:,2]**2+xp3[:,3]**2)
    (temp1,bns)=np.histogram(xp1[:,0],bins=bs,weights=tmp1)
    (temp2,bns)=np.histogram(xp2[:,0],bins=bs,weights=tmp2)
    (temp3,bns)=np.histogram(xp3[:,0],bins=bs,weights=tmp3)
    binx=0.5*(bns[1:]+bns[:-1])
    den1save[it,:]=den1
    den2save[it,:]=den2
    den3save[it,:]=den3
    avl1save[it,:]=avl1/den1
    avl2save[it,:]=avl2/den2
    avl3save[it,:]=avl3/den3
    tmp1save[it,:]=temp1/den1
    tmp2save[it,:]=temp2/den2 
    tmp3save[it,:]=temp3/den3 

CPU times: user 1min 26s, sys: 15 s, total: 1min 41s
Wall time: 4min 33s
In [9]:
%%time
plt.figure(figsize=(20,15))
plt.subplot(331);plt.imshow(den1save,origin='lower',aspect='auto',extent=[0,xmax,0,tmax],cmap='terrain_r');plt.colorbar()
plt.title(r'$density$')
plt.subplot(332);plt.imshow(avl1save,origin='lower',aspect='auto',extent=[0,xmax,0,tmax],cmap='terrain_r');plt.colorbar()
plt.title(r'$mean velocity$')
plt.subplot(333);plt.imshow(tmp1save,aspect='auto',origin='lower',extent=[0,xmax,0,tmax],cmap='terrain_r');plt.colorbar()
plt.title(r'$temperature$')
plt.subplot(334);plt.imshow(den2save,aspect='auto',origin='lower',extent=[0,xmax,0,tmax],cmap='terrain_r');plt.colorbar()
plt.subplot(335);plt.imshow(avl2save,origin='lower',aspect='auto',extent=[0,xmax,0,tmax],cmap='terrain_r');plt.colorbar()
plt.subplot(336);plt.imshow(tmp2save*25,aspect='auto',origin='lower',extent=[0,xmax,0,tmax],cmap='terrain_r');plt.colorbar()
plt.subplot(337);plt.imshow(den3save,aspect='auto',origin='lower',extent=[0,xmax,0,tmax],cmap='terrain_r');plt.colorbar()
plt.xlabel(r'$x/(c/\omega_{e})$')
plt.ylabel(r'$t\Omega_{ci}$')
plt.subplot(338);plt.imshow(avl3save,origin='lower',aspect='auto',extent=[0,xmax,0,tmax],cmap='terrain_r');plt.colorbar()
plt.subplot(339);plt.imshow(tmp3save*25,aspect='auto',origin='lower',extent=[0,xmax,0,tmax],cmap='terrain_r');plt.colorbar()
plt.show()

CPU times: user 6.71 s, sys: 140 ms, total: 6.85 s
Wall time: 6.84 s

Linear dispersion relation

In [10]:
dx=0.1
dt=dx
isav=32
nx=1024
nt=4096
kmin=2*mt.pi/(dx*nx)*(-nx/2);kmax=2*mt.pi/(dx*nx)*(nx/2)
wmin=2*mt.pi/(dt*nt*isav)*(-nt/2);wmax=2*mt.pi/(dt*nt*isav)*(nt/2)
fftbysave=fftshift(fft2(bysave[:,::-1]))  #switching the raw
#plt.pcolor(np.abs(fftexsave),cmap='terrain')    
plt.imshow(np.abs(fftbysave),origin='lower',extent=[kmin,kmax,wmin,wmax],aspect='auto',cmap='terrain')    
plt.axis([0,kmax/10,0,wmax])
plt.colorbar()
plt.clim(0,100)
plt.xlabel(r'$k(c/\omega_{pe})$')
plt.ylabel(r'$\omega/\omega_{pe}$')
plt.show()

Nonlinear evolution (mode decomposition)

In [11]:
%%time
nt=4096
nx=1024
plt.figure(figsize=(12,8))
fftbz=fftshift(fft2(bzsave[:,:]))
bzpls=np.zeros((nt,nx), dtype=np.complex)
bzmin=np.zeros((nt,nx), dtype=np.complex)
plt.subplot(131)
plt.imshow(bzsave,aspect='auto',origin='lower',extent=[0,xmax,0,tmax],cmap='bwr',clim=(-0.1,0.1))
plt.colorbar()
plt.title(r'$B_z^{raw}$')
plt.xlabel(r'$x/(c/\omega_{e})$')
plt.ylabel(r'$t\Omega_{ci}$')

plt.subplot(132)
bzpls[0:nt//2,nx//2:nx]=fftbz[0:nt//2,nx//2:nx]
bzpls[nt//2:nt,0:nx//2]=fftbz[nt//2:nt,0:nx//2]
bzpls=ifft2(ifftshift(bzpls))
plt.imshow(np.real(bzpls),aspect='auto',origin='lower',extent=[0,xmax,0,tmax],cmap='bwr',clim=(-0.1,0.1))
plt.colorbar()
plt.title(r'$B_z^+$')

plt.subplot(133)
bzmin[0:nt//2,0:nx//2]  =fftbz[0:nt//2,0:nx//2]
bzmin[nt//2:nt,nx//2:nx]=fftbz[nt//2:nt,nx//2:nx]
bzmin=ifft2(ifftshift(bzmin))
plt.imshow(np.real(bzmin),aspect='auto',origin='lower',extent=[0,xmax,0,tmax],cmap='bwr',clim=(-0.1,0.1))
plt.colorbar()
plt.title(r'$B_z^-$')

plt.show()

CPU times: user 3.76 s, sys: 551 ms, total: 4.31 s
Wall time: 4.31 s

Phase-space velocity distribution

In [12]:
it=250
iphs=128
bs=200
f1=open('xp1_%05d' %(it),'rb') 
f2=open('xp2_%05d' %(it),'rb')
f3=open('xp3_%05d' %(it),'rb')
xp1=np.fromfile(f1,dtype='float64').reshape(-1,4)
xp2=np.fromfile(f2,dtype='float64').reshape(-1,4)
xp3=np.fromfile(f3,dtype='float64').reshape(-1,4)
xion=np.r_[xp2,xp3]
plt.hist2d(xion[:,1],xion[:,2],bins=bs,cmap='terrain_r')
plt.xlabel(r'$v_x/c$')
plt.ylabel(r'$v_y/c$')
plt.title(r'$t\Omega_{ci}=%.2f$' %(it*iphs*dt/rm*wce))
plt.show()

Input parameters

In [13]:
!cat input.f90

module input
 integer(kind=4), parameter :: ng = 1024, npcel = 200, nsp =3
 integer(kind=4), parameter :: ng1 = ng+1, maxnp = ng*npcel
!
 integer(kind=4), parameter :: irun  =    1
 integer(kind=4), parameter :: nt0   =  65536*2 
 integer(kind=4), parameter :: isav  =  32
 integer(kind=4), parameter :: iphs  =  128
 integer(kind=4), parameter :: idata =  2048
!
 real(kind=8),    parameter :: al    = 102.4D+00
 real(kind=8),    parameter :: wce0  = -.2d-0
 real(kind=8),    parameter :: theta = 0.D+00
 real(kind=8),    parameter :: rm    = 25.D+0
 real(kind=8),    parameter :: dx=al/ng, dt=dx
!
 logical,         parameter :: parentwave = .false.
 integer(kind=4), parameter :: mno1  = 512
 real(kind=8),    parameter :: b01   = 1.D-00
 real(kind=8),    parameter :: v01   = 0.676581209040969D-00
 real(kind=8),    parameter :: v02   = 0.619917185536019D-01
 real(kind=8),    parameter :: freq1 = 0.198130655636275D-00
!
 real(kind=8), parameter :: wpe=1.d0, qme=-1.d0, qmi=-qme/rm
 real(kind=8), parameter :: pi=3.141592653589793d0, twopi=2*pi
 real(kind=8), parameter :: rd=theta*pi/180.d0, ak0=twopi/al
!
 contains
 subroutine com_vel(va,ve,temp)
 implicit real(kind=8)(a-h,o-z),integer(kind=4)(i-n)
 real(kind=8), dimension(nsp) :: va,ve,temp
!
 va(1) = 0.024D-00;  ve(1) = 0.00D-00;  temp(1) = 1.D-02
 va(2) =-0.00D-00;  ve(2) = 0.00D-00;  temp(2) = 1.D-03/rm
 va(3) = 0.12D-00;  ve(3) = 0.00D-00;  temp(3) = 2.5D-05
 !va(4) = 0.4242D-00;  ve(4) = 0.2242D-00;  temp(4) = 4.D-04
!
 return
 end subroutine com_vel
!
 subroutine com_den(dn,q)
 implicit real(kind=8)(a-h,o-z),integer(kind=4)(i-n)
 real(kind=8), dimension(nsp) :: dn,q
!
  dn(3) = 0.2D-00
!  dn(4) = -0.01D-00
!
 dn(1)=1.d0
 dn(2)=1.d0
 q(1)=al*wpe*wpe/(maxnp*qme)
 q(2)=al*wpe*wpe/rm/(maxnp*qmi)
 if(nsp > 2) then
   do isp=3,nsp
     if(dn(isp) < 0.d0) then
       dn(1)=dn(1)+dn(isp)
       q(isp)=q(1)
     else
       dn(2)=dn(2)-dn(isp)
       q(isp)=q(2)
     endif
     dn(isp)=abs(dn(isp))
   enddo
 endif
!
 return
 end subroutine com_den
!
end module input